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Time-dependent particle acceleration in supernova 
remnants in different environments 
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ABSTRACT 

We simulate time-dependent particle acceleration in the blast wave of a young super- 
nova remnant (SNR), using a Monte Carlo approach for the diffusion and acceleration 
of the particles, coupled to an MHD code. We calculate the distribution function of the 
cosmic rays concurrently with the hydrodynamic evolution of the SNR, and compare 
the results with those obtained using simple steady-state models. The surrounding 
medium into which the supernova remnant evolves turns out to be of great influence 
on the maximum energy to which particles are accelerated. In particular, a shock going 
through a p oc r~ 2 density profile causes acceleration to typically much higher energies 
than a shock going through a medium with a homogeneous density profile. We find 
systematic differences between steady-state analytical models and our time-dependent 
calculation in terms of spectral slope, maximum energy, and the shape of the cut-off 
of the particle spectrum at the highest energies. We also find that, provided that the 
magnetic field at the reverse shock is sufficiently strong to confine particles, cosmic 
rays can be easily re-accelerated at the reverse shock. 

Key words: MHD — ISM: supernova remnants — acceleration of particles 



1 INTRODUCTION 

Cosmic rays with energies up to at least ~ 10 15 eV are 
thought to originate in supernova remnants (SNRs). They 
are high-energy particles that have a simple power-law en- 
ergy spectrum that extends over five decades in energy. The 
need for an efficient acceleration mechanism in SNRs has 
motivated the development of the theory of diffusive shock 
acceleration (DSA), accordi ng to which particles are ac cel- 
erated at shock fronts (see iMalkov fc O'C Drury||200ll , for 
a comprehensive review). 

Young supernova remnants are ideal locations for study- 
ing this process, because of the high shock velocity, and 
the presence of a few nearby young rem nants for which de- 
tailed observation s are available (e. g. iHwang et al. 1 12004 
iBamba et al]|2005l ; lAc~ero et al.ll2009l ). The presence of high- 
energy electrons spiralling in a ~ 10 u.G magnetic field has 
been established from the emission of synchrotron radiation 
from radio wavelengths to X-rays. Both the magnetic field 
strength and the typical particle energy are inferred from 
synchrotron theory and can be used to compare with theo- 
retical predictions (I Achterberg et al.|[l994l ; IVink fe Lamind 
l2003l ; IVolk et al.ll2005l ; IVinkl2005l ). Even though synchrotron 
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emission only indicates the presence of relativistic elec- 
trons, the presence of energetic protons is suggested by 
the o bservations of TeV gamma rays (e.g. lAharonian et all 
2009), and by indications of magnetic field amplification be- 
yond that what is expected from simple shock co mpression 
(|Warren et al.1 120051 ; ICassam-Chenai et all I2008D . In addi- 
tion, there are indications that the SNR blast waves are not 
simple hydrodynamic blast waves in a single-component gas. 
They behave in a way that indicates that a significant frac- 
tion of the pre-shock energy density re sides in cosmic rays 
iDecourchelle et al.ll200ol ; lHelder et al.ll2009l ). 

Various groups work on trying to get an inte- 
gral picture of the interaction between SNR shocks, 
magnetic field s, and the associ a ted p a rticle acceler- 
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Berezhko fc Volkl l200rj ; IZirakashvili fc Aharonianl 12007 
Ferrand et al] l201dh . The difficulty is the fact that the 



process is inherently non-linear. The spectral slope of the 
particles is determined mainly by the difference between 
the plasma velocities on the two sides of the shock. This 
difference depends on the compression ratio, which in 
turn is determined by the effective equation of state of the 
gas-cosmic ray mixture: the presence of cosmic rays tends to 
soften the equation of state around the shock, which in turn 
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increases the total compression. In addition, the gradient 
in cosmic ray pressure slows down the incoming flow in the 
cosmic-ray precursor to the shock. High-energy particles 
with a larger scattering mean free path probe a larger region 
around the shock and 'feel' a higher compression ratio. For 
these reasons the spectrum flattens at the high-energy end 
in fully non-linear models. An additional nonlinearity arises 
when the magnetic fields are amplified by the streaming of 
cosmic rays. 

Cosmic rays isotropise by scattering off Alfven waves. In 
gyro-resonant scattering the scattering rate depends on the 
cosmic ray energy through the slope of the spectrum of the 
scattering waves. It is often assumed that the diffusion rate is 
described by Bohm diffusion, where the diffusion coefficient 
scales as kb oc E, with E the energy of the particle. These 
waves are self-generated by the streaming of cosmic rays 
away from th e shock, through a resonant instability (e. g. 
Skilli nelll975l) and / or through a non-resonant ins tability 
jBell fc Lucekll200ll ; lBelHl2004l ; iLuo fc Melrose! [20091 ). 

Various authors focus on the feedback of cosmic rays 
onto the hydrodynamics near the shock (iMalkovl Il99 ' 



m 

10). 



iBlasi et alj|2007l ; iPatnaude et al.ll2009l ; iFerrand et al.ll201C 

The distribution function of the particles is calculated, from 
which the cosmic ray pressure can be determined. The pres- 
sure term alters the equation of state and feeds back on 
the hydrodynamics. Alternatively, a standard power law dis- 
tribution is assumed, simplifying and speeding up the pro- 
cess, making it fast enough for application on larger scales 
jEnfilin et al.ll2007l V The disadvantage of this approach is 
that - generally speaking- the time-dependence or the in- 
fluence of a complicated magnetic field geometry can not be 
taken into account. 

Our work focuses on the time- and position depen- 
dence of cosmic ray spectra in the supernova remnant, 
where we use the test-particle approximation, neglect- 
ing feedback of the particles onto the plasma. The high- 
energy particles isotropise on both sides of the shock, 
due to scattering off waves, and are accelerated by re- 
peated crossing cycles at the shock. The acceleration of 
relativistic particles can be des cribed by a set of stochas- 
tic differential equations (S DEs) (|Achterberg fc Kriilisll 19921 ; 
iKriills fc Achte rbcre 1994]) and t his has been applied suc- 
cesfully by a number of autho r s llMarcowith fc Kirk 19991; 
van der Swaluw fc Achterberd |2004| ; iMarcowith fc Cassa 



20101 ) in various hydrodynamic codes. 

We use the adaptive mesh refi ned magneto- 



hydr o dynamics (MHD) code: AMRV AC i|Keppens et all 
120031 ; I van der Hoist fc Keppensl 120071 ) as the framework 
for our particle acceleration method. Not only do we 
calculate the acceleration/deceleration of test particles due 
to compression/expansion of the flow, we also model the 
dependence of diffusion and radiative losses, while keeping 
computational costs down with the adaptive mesh strategy. 
Our approach has the advantage that it is able to also 
tackle a more complicated circumstellar density profile 
than other models. The disadvantage is (for now) having to 
neglect the nonlinear feedback of the cosmic rays onto the 
plasma. 

We describe the different models we use to calculate the 
particle spectrum in supernova remnants in § [21 In § [3] we 
will discuss the theory of diffusive shock acceleration and 
the evolution of the supernova remnant, and derive some 



analytical estimates for the expected cosmic ray spectrum. 
The method and set-up of the simulations will be described 
in § 21 In § \S\ we will describe some test models and results 
obtained with this method. We will subsequently present the 
results for the particle spectra from the SNR models in § [6] 
and conclude with a discussion and summary in § [7] 



2 SNR MODELS 

We simulate the evolution of the SNR and concurrently cal- 
culate the cosmic ray distribution in phase-space for the 
various models. In these models we evaluate the effect of: 1) : 
the approximation that is used for the diffusion coefficient, 
ii): the adopted equation of state, Hi): the density profile of 
the background into which the supernova remnant evolves, 
and iv): the strength of the magnetic field, both in the sur- 
rounding medium, and in the ejecta. Further details of the 
latter three can be found below. We summarize the entire 
set of models used in our simulations in Table [1] 



2.1 Influence of the equation of state 

The compression ratio at the shock depends on the adiabatic 
index (specific heat ratio) of the plasma. Since the expected 
slope of the spectral energy distribution of the cosmic rays 
depends on the compression ratio, by varying the equation 
of state we can test our code and see if the results change 
according to expectation. 

Besides providing a nice test for the code, varying the 
adiabatic index is also physically relevant for systems in 
which efficient cosmic ray acceleration occurs. The adiabatic 
index of the plasma parametrises the equation of state of the 
plasma. The value of 7 = 5/3 corresponds to the case where 
the gas consists of a mono-atomic gas, and is what we use 
for simulations in which the contribution of cosmic rays is 
neglected. The value of 7 = 4/3 corresponds to the adiabatic 
index of a relativistic gas. A plasma in which a large part 
of the pressure is made up by cosmic rays has an effectively 
softer equation of state, with a value of 4/3 < 7 < 5/3. 
When cosmic rays escape from the system, they may drain 
the shock region of energy, softening the equation of state 
further. 

Approximating the hydrodynamics with a lower, but 
fixed, adiabatic index is only a first step to see how the slope 
of the cosmic ray spectrum is affected by a softer equation 
of state. In reality, the cosmic ray pressure gradient in the 
shock precursor slows down the incoming flow, which in turn 
results in a concave cosmic ray spectrum if the scattering 
mean free path increases with energy. Additionally magnetic 
field amplification may partly reverse the effect by stiffening 
the equation of state. 



2.2 The explosion environment 

We consider two situations for the environment into which 
we let the ejecta expand, i): The case of a homogeneous 
medium, such as may be applicable to Type la SNe. These 
models will be referred to as interstellar medium (ISM) mod- 
els, it): A p oc R~ 2 profile for the plasma, as will arise from 
a steady stellar wind. This situation is more applicable to 
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core collapse supernovae whose progenitors have strong stel- 
lar winds. These models will be referred to as circumstellar 
medium (CSM) models. We set the density of the homo- 
geneous ISM to a value of 2.34 x 10~ 24 g cm -3 . The wind 
parameters for the p oc BT 2 CSM are chosen such that the 
volume containing an amount of mass equal to the ejecta 
mass is equal in both cases. The ejecta mass that we use is 
3.5 Mq (Sect. [4~3)l . The radius containing an equal amount 
of mass in the ISM is P 3 .5M = (3Af/47r / 9) 1/3 ~ 3 pc. We 
employ mass loss parameters that are typical for the winds 
from red supergiants, with a value for the mass loss rate of 
M — 1.71 x 10 -5 M yr _1 and a wind velocity equal to 
^wind(10 16 cm) = 10 km s" 1 . In all cases we set the initial 
temperature of the CSM/ISM and the ejecta to T = 10 4 K. 



increase of the particle momentum every time the particl e 
cycl es from up s tream to downstream and b a ck l|Bell| 1978al) . 

iKrvmskiil rtl977T): l Axford et all Jl977tl : [Belli (|l978al lbh: 
iBlandford fc Ostrikeij (|l978MDrurvl (|l983T l derive the spec- 
trum of shock-accelerated particles, based on the energy gain 
per cycle versus the probability of being advected away and 
lost from the acceleration process. The escape probability is 
given by P csc = 4Vs/rc, where V s is the shock velocity, r the 
shock compression ratio and c is the velocity of the cosmic 
rays, equal to the speed of light. The mean time At for a rel- 
ativistic particle to complete a cycle probing the upstream 
and the downstream plasm a and the corresponding boost in 
mom entum p is given by l|Lagage fe Cesarskvl Il983l ; iDrurvl 
Il983h : 



2.3 Magnetic field 

For the magnetic field we look at the case of a constant 
magnetic field strength, with values of 3, 10, and 20 uG. We 
restrict ourselves for now to parallel shocks, where the mag- 
netic field is aligned with the shock normal, for reasons men- 
tioned in § [4] In this case the magnetic field is the same on 
both sides of the shock and the diffusion coefficient depends 
only on the particle energy in the case of Bohm diffusion. 
While the code has an MHD solver and we can follow the 
magnetic field dynamically, this only becomes interesting in 
two dimensions, something we defer until a later work. We 
solve the Euler equations and parametrise the magnetic field 
strength for use in the SDEs. 

We find (see § 16. 5p that particles diffuse far enough 
downstream to reach the reverse shock such that they are 
re-accelerated, provided that the magnetic field in the ejecta 
is strong enough to confine the particles. We therefore inves- 
tigate two different cases for the magnetic field in the ejecta: 
we either parametrise the magnetic field strength such that 
it decays with radius as BT 2 as to represent the expanding 
fossil field of the progenitor, frozen into the plasma while 
conserving magnetic flux, or we keep the magnetic field fixed 
at the same value as in the ISM/CSM. In the first case, 
the magnetic field strength for reasonable progenitor fields 
is very weak, making the gyroradius (and thus the diffusion 
length-scale) of the particles very large when Bohm diffusion 
is assumed. In the latter case, the diffusion in the ejecta pro- 
ceeds at a similar rate as in the ISM/CSM, so particles that 
encounter the reverse shock can be re-accelerated, as we will 
illustrate in Sec. 16.51 



3 THEORY 

3.1 Diffusive Shock Acceleration 

At the forw ard shock of SNRs, and possibly also at the re- 
verse shock (jHelder fc Vinkll2008l : IZirakashvili fc Aharonianl 
2010), particles are believed to be accelerated to relativistic 
energies, up to 10 15 eV for protons. The favoured mechanism 
is diffusive shock acceleration (DSA), where every time the 
particle crosses the shock interface it experiences a Lorentz 
boost, after which the momentum direction is randomised 
in the local frame. Since the plasma rest frames on either 
side of the shock differ in speed, this results in a systematic 



At = — ■ (ki + rn 2 ), 

CVs 



(1) 

(2) 



where n is the diffusion coefficient and we will use the sub- 
script 1 for upstream-, and 2 for downstream properties. 

The test particle limit results in a power-law spectrum 
in momentum, with the number density of particles per unit 
momentum following the distribution 



with 



9 = 1 + 



r + 2 



ln(l/P,. ct ) _ 

ln((p + Ap)/p) r-1 



(3) 



(4) 



the slope. Here P rct = 1 — P esc is the return probability in a 
shock crossing cycle. 

3.2 Acceleration time-scales 

The acceleration rate of relativistic protons or electrons of 
given energy E = pc^> m p c 2 ~ 1 GeV is: 

(dE\ _ ( r-l \ V 2 (5] 
V dt J dsa ~~ V 3r J (ki +r K2 ) ' 



giving a typical time-scale for acceleration time of 

= ( 1 dEy 1 _ 3r(«i +re 2 ) 
* cc -\Edt) 



(r-l)V 2 



(6) 



It will be useful to consider the case of Bohm diffusion 
for comparison with the simulations, and to look at the com- 
bined effect of acceleration and synchrotron losses, which is 
important for electrons. The change in energy at the shock 
is then given by: 

SZeBV 2 \ S B 2 E 2 

™„2 ' \'J 



(f) 



ch (1 + r)rc 



acc+syn 

for a parallel shock, where Ki = ft 2 — Db = cE/3ZeB is the 
Bohm diffusion coefficient, which depends on the gyroradius 
of the particles, B the magnetic field strength, Z the charge 
number, e the elementary electric charge and m the mass of 
the particle. 

We have introduced a term that accounts for syn- 
chrotron losses with 

\ aT -3 lo\ 

A s = oc m . (8) 

67rmc 



4 KM. Schure et al. 

Table 1. Overview of the supernova remnant models. 



MODEL: 


1SMk c 


ISMk b 


CSMk b 


ISMre B r 


CSMfc B r 


CSMk b s 


p background (g cm -3 ) 


2.34 X 1CT 24 


2.34 x 10~ 24 


oc r~ 2 


2.34 X 10~ 24 


<x r~ 2 


<x r~ 2 


diffusion coefficient 


constant 


Bohm 


Bohm 


Bohm 


Bohm 


Bohm 


magnetic field strength (u-G) 


10 


3, 10, 20 


3, 10, 20 


20 


20 


20 


reverse shock acceleration 


no 


no 


110 


yes 


yes 


no 


equation of state 


normal (7 = 5/3) 


normal 


normal 


normal 


normal 


soft (7 = 4/3) 



The mass dependence makes this term negligible for protons. 
High-energy electrons on the other hand lose a substantial 
fraction of their energy by synchrotron radiation. With a 
Thomson cross section of or = 6.65 X 1CP 25 cm 2 for elec- 
trons this gives for A s = 1.29 x 1CP 9 cm s g" 1 . The same 
loss-term can be used to account for inverse Compton losses, 
where the energy that is lost in upscattering photons from 
the microwave background corresponds to synchro tron losses 
in an equivalent magnetic field -Bcmb = 3.27 u.G (| Reynolds! 
1998). When the actual magnetic field is stronger than this 
value, synchrotron losses will dominate over inverse Comp- 
ton losses unless there is another source of photons that 
boosts the inverse Compton scattering rate. For now we ne- 
glect Compton losses, which means that we slightly over- 
estimate the maximum energy for electrons in our models 
with B ~ 3 /J.G. 

For oblique shocks the compression of the magnetic field 
and the change in residence times upstream and downstream 
have to be taken into account. 



3.3 Maximum energy of the cosmic rays 

The maximum attainable energy E max for cosmic rays de- 
pends on the size of the accelerator, the time they have spent 
there, and the strength of adiabatic- and synchrotron losses. 
In order to determine this, we follow the evolution of the su- 
pernova remnant. Initially, the ejecta expand freely into the 
surrounding medium. The expansion slowly decelerates as 
the ejecta sweep up mass from the CSM. The deceleration 
of the blast wave drives a reverse shock into the ejecta, heat- 
ing the material to millions of Kelvin and making the SNR 
prominent in X-rays. By the time the swept-up mass equals 
the ejecta mass the Sedov- Taylor phase of SNR evolution 
begins. 

Energy conservation gives us typical length- and time- 
scales, and determines the deceleration radius where the 
transition from fre e expansion phase to the Sedov- Taylo r 
phase takes place (|Sedovl Il959l ; iMcKee fc Truelovei Il995l) . 
Figure Q] shows the analytical solution for the evolution 
of the blast wave of a supernova remnant derived by 
iTruelove fe McKed [|l999T ) . We show two different cases, one 
in which the density of the ejecta is constant with radius 
(n — 0- model), and the case in which the ejecta consists of 
a constant density core with a envelope in which the den- 
sity decreases with radius as R~ 9 (n — 9-model, see also 
our description of the ejecta in the hydrodynamics in § [4j. 
The early evolution of SNR radius and blast wave velocity 
is quite different in the two cases. The e volution of the blast 
wave radius and velocity according to ITruelove fe McKeel 



(1999) for 



(9) 



ISM is given by: 

R s (t< 1st) = 1.12 P /3 

R s (i>t S T) = 1.42 (t- 0.297) 2/5 

v s (i<tsT) = 0.75 t _1/3 

v s (i>t S T) = 0.569 (1.42 (£ - 0.297)~ 3/5 ), 

with R = R/R c h and v = v/v c h- The characteristic ra- 
dius and velocity are given by f? c h = 3.07(M e j /Mg) 1 ' 3 pc, 
with M e j the ejecta mass, and v c h = R c h/tch with i c h = 
423(Mej/M Q ) 5/6 yr. The transition time is given by t S T = 
0.523. 

For a maximum residence time of the particles of i max ~ 
Rs/Vs, a simple approximation for the maximum energy for 
protons follows from Eq. [S] 



-En 



3ZeBV s R s 



(10) 



with £ CT a relation between the compression ratio of the den- 
sity and the magnetic field, where £ CT = 20 for a shock where 
the magnetic field is parallel to the shock normal, and £ CT = 8 
for a shock with the magnetic field perpendicular to the 
shock normal. The acceleration efficiency is highest around 
the Sedov- Taylor transition phase. For a simple ejecta profile 
and a homogeneous ISM (n = 0, s = 0, with n the power 
law of the density in the ejecta envelope and s the power 
law of the density of the surrounding medium) model, the 
corresponding typical radius of the blast wave and time of 
transition into this phase are given by: 



Rst 



tST 



(11) 



5/6 



-V3, 



At this point in the evolution of the SNR the maximum 
attainable energy for protons is: 



Est - 107 ( g 



-1/6 



B \ „ 



1/3 



TeV, 



(12) 



with no the number density of the surrounding medium, E51 
the explosion energy in units of 10 51 erg, and where the corre- 
spon ding radius and age of the supernova remnant are taken 
from lMcKee fc Truelovei (|l995l ). For electrons the maximum 
energy at sufficiently late times is limited by synchrotron 
losses and follows from the balance of the acceleration term 
and the synchrotron loss term in Eq. [7] This occurs when 
the electron energy equals: 



9 power law ejecta envelope into a 1 cm 



Es- 
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Figure 1. Early-time evolution of the radius (black) and velocity 
(coloured) of a SNR in a homoge neous ISM. The solid line c orre- 
sponds to the n = 9 model from iTruelove fc McKed dl999t) and 
the dashed line to the n = model. 



11.6 



V 



6000 km/s J \W^G 



B 



-1/2 



TeV, 



which corresponds to a synchrotr on loss time-scale of 
(|van der Swaluw fc Achterber3l2004l ): 



800 



B 
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VlOO TeV J 



yr. 



(14) 



In Fig. [2] we show the solution of Eq.[7]for the analytical 
estimate of the maximum energy as a function of the blast 
wave radius. From the equation we can see that the influence 
of the shock velocity on p max is very strong. The different 
evolution of the shock velocity in the n = versus the n = 9 
model therefore is reflected in the strong differences of the 
maximum energy in the two models, as can be seen in Fig. [2] 



4 METHOD 

4.1 Stochastic differential equations 

We use the method of lAchterberg fc Kriilfl (|l992f l to cal- 
culate the distribution function of cosmic rays. We describe 
the method below (details of the derivation for 2D spheri- 
cal geometry can be found in Appendix [X| and outline the 
set-up used in the hydro-simulation. 

The propagation of relativistic particles through the 
plasma can be described by a random walk. To simulate scat- 
tering off MHD waves, such as Alfven waves, the mean free 
path has a dependence on the particle momentum and the 
strength of the magnetic field. The presence of a magnetic 
field may cause the diffusion to be anisotropic, with diffu- 
sion along the field proceeding more rapidly than across the 
field. The acceleration and propagation of relativistic parti- 
cles can be described by a phase sp ace advection-diffusion 
equation (|Skillindll975l ; |joneslll99oT ): 

^ = -Vz • (UF) + Vz ■ (« • VzF) , (15) 
at 

where F(Z, t) is the particle distribution in phase space, Vz 




R (pc) 

Figure 2. Maximum energy of relativistic electrons (solid) and 
protons (dashed) for different magnetic field streng ths. The black 
lines corresponds to the n = 9 model from ITruelove fc McKeei 
(1999) and the coloured lines to the n = model. As long as the 
maximum energy is not limited by radiative losses, the influence 
of the shock velocity dominates the evolution of the maximum 
energy. The curves with the highest energies correspond to a mag- 
netic field strength of 3 u.G. The higher magnetic field strengths 
of 10, 20 and 50 uG limit the electron energy to subsequently 
lower values (normalized to the magnetic field strength). 



is the gradient operator in phase space, Z = (x, p) is the 
phase-space position vector, U = dZ/dt is the phase space 
velocity and n the diffusion tensor in phase space. 

By reordering the operators we obtain a Fokker-Planck 
equation of the standard form; 



dF(Z,t) 

dt 
Here 



= -Vz ■ [ZF(Z,t) - Vz • («F(Z,t))] 



U + Vz 



(16) 



(17) 



is an effective velocity that includes a drift term due to diffu- 
sivity gradients. Equation [16] is mathematically equivalent 
to a set of stochastic differential equations (SDEs) of the 
Ito form (|Gardinerlll99i ISaslawlll985l ; lAchterberg fc Kriillsl 
Il992h : 



dZ = Z dt + V2n ■ dW 



(18) 



The SDEs contain, apart from a regular (deterministic) term 
oc dt, a stochastic term that models the random walk. The 
Wiener process dW in the stochastic term satisfies (dWi) = 
and (dWidWj) = dijdt. The quantity *J~H = T represents 
a tensor T that should formally satisfy Ti m T mj = Kij . 

A large number of statistically independent integra- 
tions of the SDEs creates a distribution of particles in phase 
space that corresponds t o the solution o f the corresponding 
Fokker-Planck equation (|Gardinerlll99l ). 

These SDEs can be used to calculate particle ac- 
celeration in a time - depen dent flow, as first proposed in 
lAchterberg fc Kriillsl dl992|) and which has been applied 
(e.g. iKriills fc Achterberd Il994l: iMarcowith fc Kirkl [19991: 



van der Swaluw fc Achterberd 120041 : Marcowith fc Cassd 
201(1) in various hydrodynamic codes. In 2D spherical geom- 
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etry, the set of equations that needs to be solved to update 
the position of the particles in phase space (as we derive in 
Appendix [A]) is: 



du 



dR 



RdG = 



--(V • V) dt- A s BV 1 + e2 " d * ■ 



v ° + R^Te + Roe) 



(19) 



where V is the plasma velocity, u = \n(p/mc) and p = 
cos 9. The stochastic term now contains a random number 
£i drawn from a unitary normal distribution such that its 
statistical average satisfies = and = 5ij, where 

i,j stand for R and p. Note that these equations solve for 
F = R 2 F rather than for F. 

In ID slab geometry with flow velocity V(x , t) these 
equations simplify to: 

du = ia 



(20) 



., v 1 dt — \ a B 2 \/l + e 2n dt, 

dx = (v +^ S )dt+V2 K dt£ x . 

Some caution should be exercised when translating 
these equations into numerical schemes in the case that large 
magnetic field gradients are present. When Bohm diffusion 
is assumed, gradients in the magnetic field strength induce 
gradients in the diffusion coefficient. This creates the fol- 
lowing problem: In Eg. IA7l we see that the divergence of the 
diffusion tensor shows up in the effective velocity of the scat- 
tering centre of our cosmic rays. In simple schemes the strong 
gradient in the magnetic field causes the particles to under- 
sample the shock statistically, giving rise to a momentum 
spectrum of the accelerated particles that is steeper than 
can be expected in reality. A more sophisticated numerical 
scheme is needed to make sure that the spectrum reflects 
the physics, rather than mathematical artifacts (Achterberg 
& Schure, in preparation). Its implementation will be pre- 
sented in a follow-up paper. Here we will only consider par- 
allel shocks where the magnetic field strength is constant 
across the shock. 

The resulting power-law spectrum should approach the 
analytical solution of a strong shock, where the compression 
ratio depends only on the specific heat ratio 7 of the plasma, 
and the slope of the distribution depends only on r. The 
compression ratio is given by r = pi) p\ = (7 + l)/(7 — 1), 
where p denotes the plasma density. For a value of 7 = 5/3, 
r = 4, whereas for 7 = 1.1, r = 21. In terms of u = ln(p/mc) 
one has: 



(21) 



F(u) =pF(p) ocp"^" 1 ', 
with q given by Eq. [4] 

4.2 Cosmic ray injection 

As mentioned in § [TJ the injection of cosmic rays at the 
shock front is still a poorly understood process. Obser- 
vations hint at efficient inject ion, particularly for parallel 
shocks l|Berezhko fe Vol]j|2008h . but it remains problematic 



to accelerate electrons from the thermal pool to sufficiently 
high energies for them to be picked up by the DSA pro- 
cess. Since we can only model the highest-energy end of the 
particle distribution , we le ave this problem to others (see 
ISironi fc Spitkovskvl l|2009t) for relativistic shocks), and as- 
sume that the rate of particle injection is proportional to the 
mass that is swept up per unit time by the blast wave. Since 
we treat the cosmic rays in the test-particle limit, the exact 
number that is injected is not so important, but the rela- 
tive number of injected cosmic rays with time influences the 
momentum distribution and the average acceleration time 
of the cosmic rays. 

For a homogeneous ISM, the number of particles that is 
injected as the shock expands a distance dR in radius scales 



dN(p) oc R dRS(p-p ) 



(22) 



We assume mono-energetic injection with po the injection 
momentum, which for Bohm diffusion is restricted for nu- 
merical reasons, as will be explained below (Eq. I28[l . For a 
CSM that is shaped by the stellar wind of the progenitor 
the density scales with p oc R~ 2 . In this case the injection 
proceeds uniformly with respect to the blast wave radius: 



dN(p) oc dR5(p-p ) 



4-2.1 Minimum injection momentum 



(23) 



In our models, the particles are injected with relativistic 
energies. They need to be able to cross the shock in one 
scattering mean free path, otherwise they will be adiabat- 
ically accelerated and we will not see the desired effect of 
DSA. The particles are all injected with the same (relativis- 
tic) energy and the power law naturally follows provided the 
steps used for the integration of the SDEs satisfy 



V s Atdsa < Ax B < Azdiff = V2k At d 



(24) 



Here Atdsa is the time step used to integrate the SDEs, 
and Aa; a is the width of the shock transition, typically a 
few grid cells. In the case of Bohm diffusion, the particles 
are inserted with an energy that makes it possible to meet 
the criterion in Eq. [24] for the simulation of DSA from the 
earliest injection time (to = 31 yr) onwards. This puts the 
following constraint on the diffusion coefficient: 



k > 



Axj 
2Ai ri! 



with 



pc 



(cgs), 



(25) 



(26) 



3 3ZeB 

where Z is the charge number of the particle, which we take 
to be 1 since we only consider protons and electrons. The 
requirement of the advective step, Aa; a d v < Ax s , in practise 
translates to a time-step requirement 

Ax s 



Atdsa < 



(27) 



Since the shock thickness is typically resolved in about 5 
grid cells, the time step for calculating the diffusion can 
be larger than that used for the MHD time step. This one 
is limited by the courant condition where we normally use 
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Afhydro = 0.8A:Eg r id/u. We therefore 'supercycle' (i.e. in- 
crease the timestep) the diffusion of the particles to save 
computational time, and calculate the appropriate time step 
Afdsa for diffusion throughout the simulation. 

The constraint on the timestep can be used to evaluate 
the valid range for the diffusion coefficient by combining 
Ea,ns. l2"5"l to l27l The minimum injection energy for this model 
is then given by: 



3qBn 3qBv 



„At d 



2c 2 



(28) 



The scaling of the injection energy with magnetic field 
also means that our resolution requirement depends on the 
magnetic field strength assumed in the simulation. We want 
the injection energy to be at least as low as ~ f TeV, such 
that it is well below the maximum energy as can be roughly 
determined by Eq.[T3] At lower energies, the spectrum inside 
the SNR will be a featureless power law as the scattering 
mean free path of the particles is much less than, say, the 
radius of curvature of the shock, and losses are unimportant. 
The maximum number of particles injected at the end of 
the simulation is fO 6 . Particle splitting is applied when the 
particle energy increases with a factor of e, in order to reduce 
Poisson noise in the distribution at high energies. 



4.3 The hydrodynamics 

We have incorporated the equations for calculating diffu- 
sive shock acceleration into the AMRVAC framework. This 
allows us to solve the particle spectrum concurrently with 
the (magneto-)hydrodynamics of the flow, in this case the 
evolution of the supernova remnant. We model the sys- 
tem in ID spherical coordinates. We use a radial range of 
3.0 x f 19 cm (~ f pc) and a base grid of f 80 cells. Refine- 
ment of a factor two is applied dynamically where strong 
density and velocity gradients are present, up to 7 to f lev- 
els. The effective resolution is then 46080 — f84320 cells, or 
1.6xf0 14 -2.6xf0 15 cm. This is sufficient for our simulations 
of the system with B up to 20 \iG. The Euler equations are 
solved conservatively using a TVD LF scheme with minmo d 
limiting on the primitive variables dToth fc OdstrciU fi 996). 
Although this particular scheme does not resolve the shock 
within the least possible number of grid cells, the scheme is 
robust and does not have problems in dealing with the ini- 
tial conditions for the supernova ejecta, such as large density 
gradients. 

The supernova ejecta are inserted in the inner 0.1 pc 
of the grid. The density and velocity pr ofile are modelle d 
according to the self-sim ilar solution of IChevalierl (|l984h ; 
iTruelove fc McKed (|l999h . The velocity increases linearly to 
the outer edge of the ejecta, whereas the density profile con- 
sists of a constant density core and a powerlaw envelope 
with a slope of n = 9, interpolated as: 



Pej(r) 



1 + (r/r COIC ) n ' 



(29) 



The core radius and central density are such that the mass 
and energy of the ejecta are respectively M c j = 3.5 Mq and 



E Bi = 10 bl erg. 



In Fig. [3] we show the density and velocity profile for 
a SNR, 600 yr after explosion into a homogeneous medium. 
In this simple one-dimensional model four different regions 
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Figure 3. Density and velocity profile of supernova ejecta at 
a time t = 600 yr after explosion into a homogeneous medium. 
The shown profiles are for a plasma with an adiabatic index of 
7 = 5/3, yielding the expected compression ratio r = 4 at the 
shocks. 



can be identified: From inside out first we encounter the 
freely expanding ejecta, separated from the shocked ejecta 
by the reverse shock. Outside the contact discontinuity is the 
compressed ISM and the blast wave that separates the SNR 
from the undisturbed ISM. In Fig. U we show the same for 
the supernova remnant in a CSM. The evolution of the SNR 
is substantially different in the two cases. The separation 
between the blast wave and the ejecta is larger in the CSM 
model. This is due to the deceleration of the blast wave, 
which initially is high for the CSM model, but soon becomes 
slower than for the ISM background, as can also be seen in 
Fig. 1111 which will be discussed later in Sect. 16.31 



5 TEST MODELS 

In this section we will present the results from some test 
cases, in order to compare the numerical results with the 
analytical ones, and to determine the dependence of the 
spectrum on the chosen geometry. 



5.1 Analytical shock profile 

We will first present the results from a calculation where 
we describe the shock using hypertangent profile in plane 
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Figure 4. Density and velocity profile of supernova ejecta at a 
time t = 600 yr after explosion into a p oc R~ 2 medium. Com- 
pared to the evolution of a SNR in a homogeneous ISM, the re- 
verse shock is much farther inwards and the distance between the 
contact discontinuity and the forward shock much larger. 

parallel geometry. In the shock's reference frame the velocity 
V(x) along the shock normal is 

"w - (^)-(^)"( 2 Sr)- <»» 

The shock has a thickness Aa; s (in our case: 3.75 x 1CP 2 in 
arbitrary units) and is located at a location of x = x s (in 
our case at x = 9). The shock compression ratio is chosen 
to be r = 4, the value for a strong shock in a 7 = 5/3 gas. 
The time step is chosen to satisfy the condition in Eq. 1241 
At = 1.35 x 1CP 2 and the shock velocity is normalized to 
V s = 1. The diffusion coefficient is constant: n = 0.28. In 
this case, the numerically obtained cosmic ray distribution 
should closely match the analytical one. 

In Fig. [5] we show the spectrum of the particles that 
are located near the shock, and that of all the particles. The 
number of particles that is introduced at the shock increases 
linearly with time, and they are all introduced with the same 
energy. The maximum number of particles used in these sim- 
ulations is 10 6 . The time evolution shows that it takes time 
before the numerical slope approaches the analytical steady- 
state result, which in our figures would correspond to a hor- 
izontal line (meaning q — 2 for our chosen compression ratio 
of r = 4). While the spectrum at the shock indeed approx- 
imates the expected result, the spectrum of the whole par- 
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Figure 5. Spectrum as a function of time for accelerated particles 
in slab geometry and with a hypertangent velocity profile. The top 
panel shows the total spectrum, the bottom panel the spectrum 
at the location of the shock. 

ticle population is steeper. This is to be expected, because 
the overall spectrum consists of a superposition of spectra 
with different maximum acceleration times (for which Fig.JS] 
lower panel, gives a sample), and hence different cut-off en- 
ergies. The overall spectrum is therefore steeper than the 
spectrum at the shock front. 

5.2 Effects of the shock geometry: planar and 
spherical 

In this test model we explore the difference on the slope of 
the energy spectrum of the particles that arises as a result 
of the adopted geometry. We will compare this to the an- 
alytical model, which assumes a steady state solution for a 
plane parallel shock. We set up the supernova ejecta in a 
homogeneous density medium and let it evolve for 1500 yr. 
We artificially induce slab or spherical geometry, taking into 
account the volume- and surface elements depending on the 
chosen geometry, such that adiabatic losses are properly ac- 
counted for. In a SNR the radius of curvature is much larger 
than the typical diffusion path of a particle, and slab ge- 
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Figure 6. Proton spectrum for the case of a constant diffusion 
coefficient. The black line shows the spectrum for a shock in plane 
parallel geometry. The coloured line shows the case for spherical 
geometry. The maximum energy extends up to unrealistically high 
energies for which physically the particles would not be confined 
to the supernova remnant. 



ometry therefore is normally considered to be a good ap- 
proximation. However, we find that there are differences in 
the slope of the spectrum when spherical geometry is taken 
into account, and argue that geometry cannot simply be ne- 
glected. 

Figure [6] shows the differences in the energy spectrum 
that arise due to the choice of geometry. In plane-parallel 
geometry the spectrum is closer to the analytical test case 
as presented in the previous section. In spherical geometry, 
the spectrum is somewhat steeper, with a slope q w 2.15. 
This is to be expected since particles downstream of the 
shock experience adiabatic losses that are proportional to 
V/r. This reduces the mean energy gain they experience at 
the shock, which steepens the spectral slope. Additionally, 
the exact slope in a time-dependent calculation depends on 
the injection rate and its dependence on time. 




E (eV) 

Figure 7. Spectrum as a function of time for relativistic protons 
(top) and electrons (bottom) for a SNR in a constant density ISM. 
The maximum energy of the electrons is limited by synchrotron 
losses for a field of B = 10 u-G. The diffusion coefficient is a 
constant. The top curves are for a SNR age of: t « 1500 yr. 
The subsequent lower and lighter coloured curves are for times: 
t 300, 150, and 30 yr. 



6 RESULTS 

6.1 Influence of the choice of diffusion coefficient 

The particle spectrum that results from the diffusive shock 
acceleration process in supernova remnants depends on a 
number of factors as described by the equations in the pre- 
vious sections. In this section we show how the assumption 
for the diffusion coefficient leaves its imprint on the slope of 
the spectrum. 

6.1.1 Constant diffusivity 

First, we consider the case of proton and electron accelera- 
tion with a constant diffusion coefficient at a spherical shock 
that expands into a constant density ISM. We fix k to the 
value for which the particles satisfy criterion Q240 for the dif- 
fusion length. In Fig. [7] we show the integrated distribution 
of protons and electrons as function of energy, represented 



as p 2 F(p), for different times. At late times, the solution 
well below the cut-off energy should approach the analytical 
solution (corresponding to a horizontal line in this represen- 
tation). The reason why the slope remains steeper is because 
spherical geometry is taken into account, something that is 
not considered for the analytical solution. As found in the 
previous section, the spectral slope of cosmic rays in spher- 
ical geometry approximates 2.15 rather than 2. The same 
slope is found here. 

For protons, the typical acceleration time-scale now 
only depends on the shock velocity (Eq. [5}. For electrons, 
the maximum energy is limited by the balance between ac- 
celeration rate and loss rate. The acceleration rate decreases 
with the shock velocity as V 2 and, hence, with time. The loss 
rate increases with energy. The effect is therefore strongest 
at late times, which can be seen in Fig. [7] where the late- 
time curves have the maximum electron energy at a lower 
value than the early-time curves. 
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6.1.2 Bohm diffusion 

When Bohm diffusion with k oc E is assumed, the acceler- 
ation time-scale increases, for a given shock velocity, with 
particle energy as t acc oc E (Eq. [S} . ft therefore takes longer 
for the higher energy particles to assume the equilibrium 
power law slope predicted by steady-state calculations. As 
can be seen in Fig. [S] this is visible in the spectrum as a 
smooth roll-over of the spectrum beyond a certain energy. 
For sufficiently 'old' electrons synchrotron losses lower the 
maximum energy and lead to a steeper cut-off than for pro- 
tons. 

The different curves show the spectrum for different 
ages of the SNR. While for protons the cut-off energy 
steadily increases with time, this is not the case for the 
electrons. As long as the electron spectrum is determined 
by the time available for acceleration (the age of the sys- 
tem), it closely resembles the proton spectrum. However, 
once the electrons reach an energy where the synchrotron 
loss time-scale becomes comparable to the acceleration time- 
scale and/or the age of the system, the spectrum deviates 
more and more from the proton spectrum, resulting in a 
lower cut-off energy and a steeper roll-over. 

The spectral slope at the lower-energy end of the elec- 
tron and proton spectrum is the same and does not depend 
on the choice of diffusion coefficient, ft equals q ~ 2.15, the 
same value we have found in the test-case of the hypertan- 
gent shock profile in spherical geometry. This is to be ex- 
pected, since the injection rate in both the test-case and the 
supernova remnant ISM models is the same, and adiabatic 
losses operate in the same manner. 

In the case of a constant diffusion coefficient, k was 
set to the value equal to the Bohm diffusion coefficient for 
the lowest energy particles (at p = po). The acceleration 
time for the high-energy particles is therefore relatively short 
compared to the Kb models, moving the cut-off to higher 
energies. 



6.2 Choice of equation of state 

In Fig. [9] we show the energy spectrum for protons and elec- 
trons accelerated at a shock in a gas with adiabatic index 
7 = 4/3. The shock propagates into a CSM with a p oc R~ 2 
density profile. We represent the spectrum as p 3 ^ 2 F(p) and 
compare the results to that of a 7 = 5/3 plasma model. 
For 7 = 4/3 the expected spectral slope at a steady shock is 
q — 1.5 (see Eq.|4]). We see that the simulations bear out this 
expectation. The other model parameters are B — 20 |J.G, 
imax ~ 1500 yr, and a CSM background. For a shock prop- 
agating into a constant-density ISM the same change in the 
spectral slope will result. 

Varying 7 is a nice test to see if the slope changes ac- 
cording to expectations, but also is a way of determining in 
which direction results change when cosmic ray acceleration 
is efficient enough to change the equation of state. In real- 
ity, the cosmic ray pressure that causes the softening of the 
equation of state is localized in the region around the shock, 
where most of the relativistic particles are located. We do 
not include this position dependence and keep the adiabatic 
index fixed to the same value over the entire grid. 

In our approach the softer equation of state results in a 
higher compression at the shock and a slower blast wave, as 





Figure 8. Spectrum as a function of time for relativistic protons 
(top) and electrons (bottom) for a SNR in a ISM. The maximum 
energy of the electrons is limited by synchrotron losses for a field 
of B = 10 |J.G. The approximation of Bohm diffusion is used for 
the scattering of the particles. As can be seen from the abrupt 
cut-off at the low-energy end of the spectrum, the injection energy 
of the particles is around 1 TeV. 



can be seen in Fig. 1101 The cut-off energy for protons there- 
fore lies at a lower value than in the 7 = 5/3 model. Since 
we assume a fairly high magnetic field in this simulation and 
the 7 = 5/3 case we use for comparison, the electron energy 
is limited by synchrotron losses. Therefore in both models 
electrons have a similar value for E nlax . 



6.3 CSM versus ISM models 

In this section we compare the results from the CSM and 
ISM models with B = 3 \iG. Bohm diffusion is assumed 
in both cases. As we already saw in the analytical calcu- 
lations in Sect. 13.21 the maximum particle energy depends 
strongly on the shock velocity V a . The density of the medium 
in which the SNR expands affects the shock velocity and 
therefore leads to differences in the cosmic-ray acceleration 
rate between the CSM and the ISM models. 

A blast wave expanding into the CSM hits a relatively 
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Figure 9. The proton (solid line) and electron (dashed line) 
spectra as resulting from a SNR evolving into a CSM with a 
softer equation of state: the CSMkb20s model (coloured lines). 
The maximum energy of the electrons is limited by synchrotron 
losses because of the high magnetic field strength (B = 20 |xG). 
The black lines indicate the spectra for the same model but with 
7 = 5/3 
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Figure 10. The radius of the blast wave (solid line) and the 
reverse shock (dashed line) is plotted for the CSM model with 
7 = 4/3 (coloured lines) and compared to the CSM 7 = 5/3 model 
(black). The higher compression ratio of the 7 = 4/3 model leads 
to the shorter distance between the blast wave and the reverse 
shock. The blast wave radius is typically smaller and thus its 
velocity, too, resulting in a lower -Emax for the protons accelerated 
in this simulation. 



dense medium early on in its evolution. As a result, the 
initial velocity is smaller but the deceleration proceeds more 
slowly, as the swept-up mass increases with radius as M sw oc 
R, as opposed to M sw oc _R 3 in the ISM case. Ultimately this 
results in a shock with a higher velocity at the end of the 
simulation (1500 yr). In the ISM model the initial shock 
velocity is higher, but the deceleration much more severe, 
and the maximum attainable particle energy is lower, at 
least for the model parameters used here. In Fig. [TT]we show 
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Figure 11. Top: evolution of the radius of the forward (solid line) 
and the reverse (dashed line) shock of the SNR in the CSM (black) 
and ISM (yellow) models. Bottom: the corresponding velocity of 
the blast wave and the reverse shock (relative to the upstream 
velocity). 



the evolution of the radius of the forward and the reverse 
shock in the top panel, and in the bottom panel the velocity 
of the blast wave for the SNR in CSM, ISM. 

The injection rate of particles is taken to be propor- 
tional to the mass that is swept up per unit time by the 
blast wave. As a result the age distribution of cosmic rays in 
the CSM and ISM models differs. This difference affects both 
the maximum energy and the shape of the overall spectrum. 
This is shown in Fig. 1121 The proton/electron spectrum in 
the CSM model, represented as p 2 F(p), is slightly concave. 
This arises because of the higher fraction of 'old' particles in 
the CSM model compared to the ISM model. On average, 
these older particles have a higher energy and (for Bohm 
diffusion) a larger diffusivity. Low-energy particles on the 
other hand are more rapidly swept away from the shock. 
This simulation shows the differences that arise in a time- 
dependent calculation with respect to the results at a steady 
(unchanging) shock. 

In Fig. [13] we show the maximum energy of the par- 
ticles as extracted from the simulations. Since the slope of 
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Figure 12. Spectra for protons (solid) and electrons (dashed) for 
CSM (black) and ISM (coloured) models with B = 10 ^G. 



the overall spectrum in this case is about q — 2.15, we de- 
fine -E m ax as the e-folding energy, where p 2 ' 15 F(p) for the 
cumulative spectrum decreases to a value smaller than 1/e 
times the value at lower energies. The higher average veloc- 
ity of the blast wave when it evolves into a CSM increases 
Emax by a factor 2 — 4 for the CSM models. Due to the low 
magnetic field strength, the synchrotron loss time for this 
model is significantly longer than the running time of the 
simulations (~4x 10 4 yr versus 1500 yr). Therefore there is 
no significant difference between the proton and the electron 
spectrum. In figure [13] we overplotted the maximum energy 
as calculated analytically in Sect. U for protons (since, be- 
cause of the low magnetic field strength in this particular 
model that for electrons is about the same). 

There is a clear difference between -E max derived from 
the analytical estimate with that from the simulations. For 
the CSM, the analytical model underestimates this value, 
whereas for the ISM the analytical estimate is significantly 
larger than the value we derive from the simulations. This 
conclusion holds when we calculate -E max as the e-folding 
energy for the spectrum at the shock, but also for the spec- 
trum of all particles in the SNR. We attribute this to the 
different age distributions in the CSM and ISM models, as 
explained above. 



6.4 Results for different magnetic field strengths 

In the previous section we showed results for a magnetic field 
strength of 3 [iG. While this value is too low to induce sig- 
nificant differences between proton and electron spectra in 
these models, differences arise for stronger magnetic fields. 
The maximum proton energy increases as E m&K oc B in the 
radiation loss-free case that applies here. The electron en- 
ergy decreases roughly as -E ma x oc 1/y/B for electrons if 
synchrotron losses are important. 

In Fig. [14] we compare the maximum particle energy for 
three different magnetic field strengths (3, 10, and 20 \iG) for 
the case of a CSM (dashed curves) and the ISM model (solid 
curves) . Bohm diffusion is assumed in both cases. For clarity, 
we plot the energy divided by the magnetic field strength to 
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Figure 13. -E max for protons (solid) and electrons (dashed) for 
CSM (black) and ISM (coloured) models with k b and B = 3 \iG. 
The analytical solution for the CSM (ISM) model (without ra- 
diation losses) is plotted with dash-dot-dots in black (coloured). 



scale away the loss- free -E ma x oc B behaviour. From Eqs. [13] 
and [TJ] we find typical synchrotron loss time-scales of 7 x 
10 3 yr for B = 10 \iG, and 2 x 10 3 yr for B = 20 uG for 
the ISM models. For the ISM model with B = 20 uG we see 
a significant difference between the maximum proton and 
electron energy. For the ISM models with lower magnetic 
field strengths, synchrotron losses are not very important 
for the energies reached in the simulations. For the CSM 
models, however, the maximum cosmic ray energy at a given 
shock radius is higher due to the larger shock velocity. The 
higher energies decrease the synchrotron loss time-scales for 
electrons. As can be seen in Fig. 1141 (both for the B — 10 |J.G 
and the B — 20 [iG case) the maximum electron energy 
is essentially determined by synchrotron losses for a shock 
radius R > 4 pc. 

When electron acceleration operates in a regime where 
the synchrotron losses roughly balance the acceleration, 
we can compare -E max to the exact asymptotic solutions 
for the cut-off region of the spectrum in the case of a 
steady, plane parallel shock, as derived analytically by 
IZirakashvili fc Aharonianl (|2007l ). In our notation: 



-E ma x — 



9ire m c cV s 



a T B (<7 + 2)' 



(31) 



At early times, when -B max is age-limited, this solution is 
not valid. Later (typically for a shock radius R > 3 — 4 pc) 
we find this steady-state result slightly underestimates -B max 
for electrons. Asymptotically, our value for the electron -E max 
seems to converge to the steady state result. 

We illustrate the difference between age-limited and 
loss-limited cosmic ray acceleration in Fig. 1151 There we 
show the cosmic-ray distribution of in energy-radius phase- 
space for a SNR with an age of ~ 1300 yr, using the CSM 
model with B = 20 iiG. The relativistic particles are accel- 
erated at the forward shock, where the particles reach the 
highest energy. The maximum energy of protons is limited 
by the age of the remnant, and therefore still increases at 
the end of the simulation. 
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Figure 14. Maximum energy as a function of blast wave radius 
for relativistic electrons (dashed) and protons (solid). The top 
(bottom) panel shows the case for the SNR in a ISM (CSM), 
for values of the magnetic field of B = 3, 10, 20 uG, assuming 
Bohm diffusion . The dotted lines show the stea dy-state solution 
for electrons by[Zirakashvil i~fc Aharonianl (2007) for the different 
magnetic field strengths. 



Figure 15. Cosmic ray protons (top) and electrons (bottom) for 
the CSMk s 20 model for a SNR age of ~ 1300 yr. The diagonally 
shaded area indicates the region > 4.6 diffusion lengths upstream 
of the shock, as calculated with the advection time-scale. The 
horizontally shaded area shows the region limited by 4.6 diffusion 
lengths when the time-scale is determined by synchrotron losses. 
The dotted line shows the evolution of p max as a function of radius 
and the thick black dot indicates p ma x at the shock. 



Sufficiently far downstream from the shock, the 'older' 
cosmic rays are located that have been advected away from 
the shock. The maximum energy of this older population is 
lower, as reflected in the distribution of particles. The dotted 
line in the figure shows the cut-off energy as a function of 
radius, and the thick black dot marks the intersection of the 
shock location with the cut-off energy in the energy-radius 
plane. The majority of cosmic rays is located between the 
forward shock and the contact discontinuity. 

Upstream, cosmic rays diffuse ahead of the shock over a 
typical distance set by the diffusion length Ld ~ V s /kb(E) oc 
E. The diagonally shaded area is the region more than 4.6 
diffusion lengths ahead of the shock, in which 99 per cent of 
the particles should be located (F(p) = F(p) f e - Ax/x ° = 
0.01F(p)o). In this model the attainable electron energy is 
limited mostly by synchrotron losses, and considerably less 
than the proton energy at an age of ~ 1300 yr. The high- 
energy electrons therefore diffuse over smaller distances com- 



pared to protons. The horizontally shaded region in the elec- 
tron plot shows the region more than 4.6 loss-limited diffu- 
sion lengths ahead of the shock, (Z/i OS s,d = \/2nti oss , with 
ti oss = Qirmlc 3 /a T B 2 E). 

In Fig. [16] we follow the population of protons that 
were injected during the first 300 yr of the simulation, 
and track their distribution in the energy-radius plane as 
the simulation progresses. At the end of the simulation (at 
t « 1500 yr), they constitute about 10 per cent of the total 
proton population. The maximum energy in the distribution 
of these protons turns out to be similar to that of the entire 
proton poulation. This is due to the fact that for this partic- 
ular model (ISMfcslOp), the Sedov- Taylor time is roughly 
equal to the time at which we start tracking the population 
of particles. As we discussed earlier, the acceleration effi- 
ciency is highest around the Sedov- Taylor transition, after 
which the maximum energy of the particles only increases 
by a factor of a few. 
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Figure 16. The evolution of the proton distribution in the 
energy-radius plane as a function of time for the ISMk^IO model. 
The filled contour plots the total particle distribution for times 
(top to bottom) 381, 635, 953, and 1587 yr. The contour lines 
track the population of particles that is injected up to time 381 
yr, after which the particles from this early population are still 



We have determined that the overall spectrum is some- 
what steeper than the canonical q — 2 power law. We at- 
tribute this to the losses that are not taken into account in 
the analytical calculations for planar, steady shocks. It is 
therefore interesting to look at the spectrum at the shock, 
one diffusion length Ld(-E max ) upstream and downstream, 
and compare it to the spectrum of all particles in the SNR. 
In Fig. [17] we show these spectra for the CSMkbIO model for 
electrons and protons. It becomes apparent that the spec- 
trum at the shock follows the analytically predicted power 
law slope quite closely. Upstream of the shock only the most 
energetic particles can be found. Since we probe the re- 
gion upstream at a diffusion length for -E ma x, this is where 
the peak of the upstream spectrum is found. At the same 
distance downstream, the slope of the spectrum is slightly 
steeper than at the shock. The spectrum is loss-limited for 
the electrons, as is evident from the sharper cut-off. 

In (semi-) analytical models, the spectrum at the shock 
in the presence of losse s is often described by (e.g. 
iMalkov fc O'C DrurvfeoOll ): 



F(p) oc p 9 e 



-(p/pn 



(32) 



with a = 1 for protons and a — 2 for electrons. We find that 
the simulated spectra quite nicely follow this exponential 
cut-off prescription for the cumulative particle spectrum. At 
the shock, the cut-off for protons is slightly sharper than 
usually assumed: we find a s=s 1.2 — 1.3. For electrons, the 
cut-off is less sharp, but still sharper than for protons: it 
closely follows a ~ 1.7. This will depend on to what extent 
the electron spectrum is terminated by synchrotron losses. 

6.5 Re-acceleration at the reverse shock 

At the forward shock, the magnetic field is likely amplified 
by the Bell-Lucek mechanism. Whether a seed magnetic field 
with values as low as to be expected at the reverse shock 
from an expanding SNR would be sufficient to trigger mag- 
netic field amplification, and whether this amplification ulti- 
mately leads to a field strength that is high enough to confine 
particles at the rev erse shock, is still a matter of discussion 
l|Ellison et alj|2005h . Some studies require higher seed values 
for the field strength, while other studies obse rve magnetic 
field amplification from an almost zero field (|Chang et al.l 
2008). We explore what happens when the magnetic field 
in the ejecta is strong enough to confine cosmic rays when 
Bohm diffusion is assumed. For a stronger magnetic field, the 
particles can undergo many shock crossings because of their 
smaller mean free path, and effectively can be re-accelerated 
at the reverse shock. A small mean free path also makes 
the expansion losses in the flow relatively less important. 
For simpliciy, we consider the case where the magnetic field 
strength in the ejecta is equal to that in the CSM/ISM. 

In our simulations we observe that very energetic cos- 
mic rays diffuse far enough downstream into the remnant to 
reach the reverse shock. Re-acceleration can then occur if 
the local magnetic field is sufficiently strong. Fig. [18] shows 
that the distribution in the energy-radius plane now has two 
peaks in energy. The highest energy (around 10 14 eV) is still 
attained at the forward shock, but acceleration at the reverse 
shock now causes an additional peak in the energy spectrum 
at its location (R ~ 3 pc for t ~ 700 yr) . 

If magnetic fields turn out to be indeed high enough to 
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Figure 17. Proton (top) and electron (bottom) spectra 
(CSMksIO model) at the shock and at locations of one diffu- 
sion length (for an energy -E max ) upstream and downstream of 
the shock, with Ax^g = y/2nt^ v and r, a( j v = re/V^. The thin 
solid line shows the cumulative spectrum. 



confine cosmic rays at the reverse shock, the picture sketched 
here is over-simplified. Apart from the cosmic rays that have 
diffused far into the SNR from the forward shock, an addi- 
tional cosmic ray component may arise if higher-Z elements 
from the supernova ejecta are accelerated at the reverse 
shock. Figure [19] shows the influence of re-acceleration on 
the shape of the spectrum. The proton spectrum is slightly 
flatter, and extends to a slightly higher energy. The electron 
spectrum is hardly modified as synchrotron losses lead to a 
cut-off around 10 13 5 eV. 

Figure [20l shows that when re-acceleration takes place, 
the maximum proton energy increases compared to the situ- 
ation without re-acceleration at the reverse shock. For elec- 
trons this initially is also the case, as long as the particle 
energy is limited by the time spent in the source. However, 
in the loss-limited regime, the maximum electron energy is 
less compared with the case of weak fields near the reverse 
shock as the larger volume with a high magnetic field leads 
to more synchrotron losses that are apparently not compen- 
sated by re-acceleration at the reverse shock. 





Figure 18. The distribution of protons (right) and electrons 
(left) for a SNR in a ISM (top) or CSM (bottom), when the 
magnetic field at the reverse shock is amplified to the same levels 
as at the forward shock. We assume a magnetic field strength of 
20 u.G, and the distribution is plotted for a remnant of an age 
of about 700 yr. The location of the shock is indicated at the 
location of -E max with a black dot. 
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Figure 19. Spectra when the magnetic field at the reverse shock 
is equally strong as at the forward shock for the CSMftg20r model 
(coloured). The solid (dashed) lines show the proton (electron) 
spectrum for a remnant age of t = 1587 yr. The black curves 
show the spectra when there is no re-acceleration at the reverse 
shock. 



7 DISCUSSION AND CONCLUSIONS 

In this paper, we have calculated diffusive shock accel- 
eration through the first-order Fermi mechanism, using 
stochastic differential equations. We treat the cosmic rays as 
test-particles and follow their acceleration and propagation 
along with the evolution of a superno va remnant. We have 
extend ed the model as desc ribed by Achterbcr g fc Kruiisl 
lll992l) (and used by e.g. Ivan der Swaluw fc Achterberd 
12004 ) to account for spherical geometry, something that is 
relevant when using this model to simulate cosmic ray accel- 
eration in supernova remnants. The model is set up generi- 
cally so it can use local magnetic field strengths to calculate 
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Figure 20. Maximum energy as a function of time for both rel- 
ativistic electrons (dashed) and protons (solid) in a SNR where 
the magnetic field at the reverse shock is 20 u.G, as it is in the 
rest of the remnant (coloured) , compared to the model with neg- 
ligible magnetic field in the ejecta (black). As in the case of no 
re-acceleration, the energy the particles obtain in the CSM mod- 
els (top) is higher than in the ISM models (bottom). While for 
protons the energy is higher in the 'reacceleration' models, for 
the electrons this turns out mostly not to be the case. Since the 
added magnetic field induces synchrotron losses also downstream 
of the contact discontinuity, the maximum energy in these models 
in fact turns out to be lower when the SNR is in its loss-limited 
regime. 



the diffusion coefficient of the test-particles. However, in this 
paper we employ a constant magnetic field strength in the 
ISM/CSM. The ejecta magnetic field is parametrised in such 
a way that it either decays as expected for an expanding field 
that is frozen into the plasma, or we fix the magnetic field 
at the same strength as used for the ISM/CSM. This is be- 
cause the numerical method that we use is not applicable in 
cases where strong gradients of the diffusion coefficient that 
may result from gradients in the magnetic field are present. 
The unique set-up of the calculation of the acceleration of 
the cosmic rays concurrent with the evolution of the SNR 
allows us to model the time- and location-dependent spec- 
trum more accurately than other models. The disadvantage 



of our method is that it does not (yet) include feedback of 
the cosmic ray pressure onto the local plasma. Our calcula- 
tions show the following results. 



Energy spectrum and spectral slope: 

With our method we can accurately model the spectrum of 
particles, where the slope is close to q = 2 for an adiabatic 
index of 5/3, and to q = 1.5 for an adiabatic index of 4/3. 
The slope of the spectrum depends on the location: near 
the shock the slope of the spectrum is closest to the ana- 
lytically predicted value for steady and planar shocks, while 
the slope of the cumulative spectrum (all particles in the 
source) is steeper. The cumulative spectrum is steeper than 
that at the shock because it also contains particle popula- 
tions from downstream of the shock, for which the particles 
were in the acceleration process for a shorter time and hence 
have lower cut-off energies. We also find that in spherical 
geometry the overall spectrum is slightly steeper when com- 
pared to a simulation in slab geometry. The main reason 
for this is the inclusion of adiabatic expansion, which causes 
the acceleration process to be slightly less effective: parti- 
cles lose a fraction of their energy in the downstream part 
of the shock crossing cycle, reducing the mean energy gain 
per cycle. The analytical solution for the steady state par- 
allel geometry where adiabatic losses are excluded therefore 
does not strictly apply in spherical geometry. The detailed 
shape of the spectrum also depends on the injection rate of 
particles as a function of time, which differs in the ISM and 
the CSM models that we employ. 



Maximum particle energy: 

There are additional differences between the results of our 
approach and those obtained using steady-state analytical 
models. The cut-off energy i5 ma x is different from that ob- 
tained analytically from the balance between acceleration 
and losses. For protons the S max is determined by the age 
of the source. The CSM simulations show a higher -E ma x 
than the analytical estimate, whereas for the ISM models 
the trend is the other way around. We attribute this to the 
difference in cosmic ray age distribution, where, since we 
assume the injection rate to be proportional to the amount 
of swept-up mass, the fraction of 'older' particles, which are 
accelerated for a longer time, is relatively high in the CSM 
model and relatively low for the ISM model. For Bohm dif- 
fusion, the high-energy end of the spectrum shows a distinct 
cut-off, either caused by the finite age of the SNR or, in the 
case of electrons, by synchrotron losses. 



Shape of the cut-off in the energy spectrum: 

The shape of the cut-off region for the cumulative spectrum 
follows quite nicely the quasi-exponential drop that is often 
assumed. However, if one looks at the spectrum in the close 
vicinity of the shock, the proton spectrum falls off slightly 
sharper than for the overall proton spectrum, while for elec- 
trons in the loss-limited regime the cut-off of the spectrum 
at the shock is more gradual than the cut-off in the overall 
spectrum. 
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CSM versus ISM: 

The strong dependence of the acceleration rate on the shock 
velocity causes the cosmic ray distribution to become sensi- 
tive to the surrounding medium. The density profile of the 
environment into which the supernova explodes determines 
the shock velocity and its evolution. The shock velocity (to- 
gether with the magnetic field) therefore determines the ac- 
celeration rate and the maximum particle energy that can 
be attained. Because in our models the average shock veloc- 
ity is much higher for a SNR expanding into a CSM, and 
because the average cosmic ray age in the remnant is larger 
in the CSM case, -E max is larger in our CSM models than in 
the ISM models. 

This suggests that core-collapse SNe in dense environ- 
ments, such as expected around a red supergiant, may be the 
most efficient particle accelerators and therefore the domi- 
nant contributors to cosmic rays up to the "knee" -energy. 
In the absence of a significant magnetic field in the ejecta, 
the particles are mostly located between the blast wave and 
the contact discontinuity. The distance between those is also 
sensitive to the velocity-evolution of the SNR and therefore 
the surrounding medium, and determines how long electrons 
are subjected to synchrotron losses. 

Re- acceleration at the reverse shock: 

If the magnetic field is sufficiently amplified at the re- 
verse shock, cosmic rays that are advected away from the 
blast wave can be re-accelerated at the reverse shock. This 
has important consequences for the maximum energy and 
the distribution of the cosmic rays. The maximum at- 
tainable energy for protons becomes significantly higher if 
re-acceleration at the reverse shock takes place, whereas 
for electrons the additional synchrotron losses in the now 
strongly magnetised SNR interior can have the opposite ef- 
fect. In reality it is conceivable that due to localized mag- 
netic field amplification, the net effect is a higher maximum 
energy for electrons, too. 

Overall, we conclude that a time-dependent calculation 
of diffusive shock acceleration in SNRs shows significant dif- 
ferences compared with steady-state plane-parallel analyti- 
cal models. The environment of the SNR has a large impact 
on the maximum attainable energy of the cosmic rays. The 
age distribution of the cosmic rays determines whether a 
time-dependent approach yields higher or lower maximum 
attainable energies. 



APPENDIX A: 
GEOMETRY 



SDES IN SPHERICAL 



In the equation for the distribution function of relativistic 
particles we need to take account of geometrical effects in 
order to describe adiabatic losses in a spherical shock ge- 
ometry. The advection-diffusion equation (Eq. I15|l can be 
written, in spherical coordinates (R , 9 , <f)) with axial sym- 
metry {d/d<f> — 0), as: 



^F(R ,9,p,t) = -Sr - S e - S p , (Al) 



with 

S R = ±^ 
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We assume isotropic spatial diffusion: k = kI, with I the unit 
matrix in configuration space. The magnetic field orientation 
is not taken into account for the diffusion rate, and diffusion 
in momentum space (Fermi-II acceleration) is neglected. 

In order to be able to apply the Ito method (Eq. I18p . 
we have to re-order the differential operators in Sr and Se, 
such that they conform to the Fokker-Planck standard form 
(cf. Eq.QD). We substitute F = R 2 F into Eq.rXT]to get: 



dF 

-q£ — —Sr — Se — S p . 



(A6) 



The operator of Sr = R 2 Sr can be rewritten to the stan- 
dard form as follows: 



d_ 

OR 

d_ 

OR 



(A7) 



OR 



This now conforms with the standard form with an effective 
radial velocity 
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This is the velocity that enters the equivalent SDE. 

In a similar fashion, Se = R 2 Se can be rewritten by 
changing the independent variable to /i = cos# (so that 
d/dO = -sindd/dfj,): 
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For we substitute u = ln(p/mc) and dp = p du to obtain 
the equivalent equation for u: 





pS p 



du 



f-V-V + A^Vl + e 2 ") F 



(All) 



We have added on the right-hand side the effect of syn- 
chrotron losses, with A s given by Eq. [8] As stated before: 
synchrotron losses may be neglected for protons but are im- 
portant for electrons. 

We now have the advection-diffusion equation in the 
standard form of the Fokker-Planck equation to which we 
can apply the Ito method to calculate the particle distri- 
bution function in a stochastic manner (Eg. I18[). We solve 
the following equations numerically to update the position 
of the particles in phase space: 

du = -~(V • V) dt - \ s B 2 yJl + e 2 " dt , 

dR = Vf dt + V2ndt £ R 

= (v R + ±?^\dt + V2^dttH, (A12) 



dfi 



Vfdt+ y/2D„dt 



Ve sin 9 dD u \ , tt; — r 

+ ) dt+ y/2D„dt^. 



R du 

If we revert from [i to 6 the last equation becomes: 
Rd6 = (Vg + 



2k 



lUtm9 + ^)*t-V^dt^ 



(A13) 



The stochastic terms should satisfy — and = 
Sij, where i,j stand for R and jj,. Note that these equations 
solve for F = R 2 F rather than for F. In a slab geometry 
with flow velocity V(x , t) along the shock normal the corre- 
sponding equations simplify to Eqn (120 ^ of the main paper. 
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